Notebook 5: Rational Method and Travel Times¶
In this notebook, we look at two ways of estimating a runoff hydrograph from precipitation data.
First, we’ll use the rational method to approximate peak flow and the maximum water level at the basin outlet, and then we’ll use an open-source library to make a higher resolution estimate of flow accumulation paths and stream networks using a digital elevation model.
# import required packages
import pandas as pd
import numpy as np
import math
# advanced statistics library
from scipy import stats
import matplotlib.pyplot as plt
import matplotlib.patches as mp
import matplotlib.colors as colors
# SEE COMMENTS ABOUT PYSHEDS LIBRARY IN NEXT CELL
from pysheds.grid import Grid
import warnings
warnings.filterwarnings('ignore')
from bokeh.plotting import figure, show
from bokeh.io import output_notebook
from bokeh.models import LinearColorMapper, LogTicker, ColorBar, BasicTickFormatter
from bokeh.io import output_notebook
output_notebook()
%matplotlib inline
Import Precipitation Data¶
For this exercise, we will use historical climate data from the Meteorological Service of Canada (MSC) station at Whistler, BC.
# calibration data
df = pd.read_csv('../../data/Whistler_348_climate.csv',
index_col='Date/Time', parse_dates=True)
# note that the 'head' command shows the first five rows of data,
# but in this case the columns are abbreviated.
# print(df.head())
# list all the columns
# print('')
# print('__________')
for c in df.columns:
print(c)
stn_name = df['Station Name'].values[0]
Station Name
Year
Month
Day
Data Quality
Max Temp (°C)
Max Temp Flag
Min Temp (°C)
Min Temp Flag
Mean Temp (°C)
Mean Temp Flag
Heat Deg Days (°C)
Heat Deg Days Flag
Cool Deg Days (°C)
Cool Deg Days Flag
Total Rain (mm)
Total Rain Flag
Total Snow (cm)
Total Snow Flag
Total Precip (mm)
Total Precip Flag
Snow on Grnd (cm)
Snow on Grnd Flag
Dir of Max Gust (10s deg)
Dir of Max Gust Flag
Spd of Max Gust (km/h)
Spd of Max Gust Flag
Plot the Data¶
It’s always a good idea to begin by visualizing the data we’re working with.
# plot flow at Stave vs. precip at the closest climate stations
p = figure(width=900, height=400, x_axis_type='datetime')
p.line(df.index, df['Total Precip (mm)'], alpha=0.8,
legend_label="Total Precip [mm]", color='dodgerblue')
p.line(df.index, df['Total Snow (cm)'], alpha=0.8,
legend_label="Total Snow [cm]", color="firebrick")
p.line(df.index, df['Total Rain (mm)'], alpha=0.8,
legend_label="Total Rain [mm]", color='green')
p.legend.location = 'top_left'
p.legend.click_policy = 'hide'
p.xaxis.axis_label = 'Date'#
p.yaxis.axis_label = 'Daily Rain [mm] / Snow[cm] Volume'
show(p)
Simplified Version of Rainfall-Runoff¶
First, isolate a single precipitation event to use for estimating a runoff hydrograph. Let’s find a nice week for skiing:
fig, ax = plt.subplots(1, 1, figsize=(16,4))
sample_start = pd.to_datetime('2014-12-01')
sample_end = pd.to_datetime('2014-12-15')
sample_df = df[(df.index > sample_start) & (df.index < sample_end)][['Total Precip (mm)', 'Total Snow (cm)', 'Total Rain (mm)']]
# print(sample_df.head())
ax.plot(sample_df.index, sample_df['Total Precip (mm)'], label="Total Precip [mm]")
ax.plot(sample_df.index, sample_df['Total Snow (cm)'], label="Total Snow [cm]")
ax.plot(sample_df.index, sample_df['Total Rain (mm)'], label="Total Rain [mm]")
ax.set_xlabel('Date')
ax.set_ylabel('Precipitation')
ax.set_title('{}'.format(stn_name))
plt.legend()
<matplotlib.legend.Legend at 0x7f1aca63a520>
First, imagine we are some unfortunate parking lot attendant working a shift in Whistler Village at Parking Lot 5, and we are told by our cruel supervisor we have stand at the lowest point of the parking lot: a catchment basin with an area of \(1 km^2\) where water runs off into FitzSimmons Creek. The sky looks angry, but we’re running late for work and put on our running shoes instead of our sturdy waterproof boots.
Next, assume the travel time is effectively zero across our entire basin (precipitation takes no time to travel to the outlet once it falls on the parking lot surface). Is this a reasonable assumption in general?
Under these assumptions, lets reconstruct a runoff hydrograph at the outlet. First, look at the precipitation data over the twelve days of the big storm.
print(sample_df)
Total Precip (mm) Total Snow (cm) Total Rain (mm)
Date/Time
2014-12-02 0.0 0.0 0.0
2014-12-03 0.0 0.0 0.0
2014-12-04 4.4 6.0 0.0
2014-12-05 12.0 0.0 6.0
2014-12-06 12.0 0.0 6.0
2014-12-07 2.8 0.0 2.8
2014-12-08 30.5 0.0 30.5
2014-12-09 89.6 0.0 89.6
2014-12-10 74.4 0.0 74.4
2014-12-11 19.2 0.0 19.2
2014-12-12 3.6 0.0 1.8
2014-12-13 0.0 0.0 0.0
2014-12-14 0.0 0.0 0.0
Convert Volume to volmeteric flow units¶
Runoff is typically measured in \(\frac{m^3}{s}\), so convert \(\frac{mm}{day}\) precipitation to \(\frac{m^3}{s}\) runoff.
# convert to runoff volume
drainage_area = 1 # km^2
# runoff is typically measured in m^3/s (cms for short -- cubic metres per second),
# so express the runoff in cms
sample_df['runoff_cms'] = sample_df['Total Rain (mm)'] / 86.4
print(sample_df)
Total Precip (mm) Total Snow (cm) Total Rain (mm) runoff_cms
Date/Time
2014-12-02 0.0 0.0 0.0 0.000000
2014-12-03 0.0 0.0 0.0 0.000000
2014-12-04 4.4 6.0 0.0 0.000000
2014-12-05 12.0 0.0 6.0 0.069444
2014-12-06 12.0 0.0 6.0 0.069444
2014-12-07 2.8 0.0 2.8 0.032407
2014-12-08 30.5 0.0 30.5 0.353009
2014-12-09 89.6 0.0 89.6 1.037037
2014-12-10 74.4 0.0 74.4 0.861111
2014-12-11 19.2 0.0 19.2 0.222222
2014-12-12 3.6 0.0 1.8 0.020833
2014-12-13 0.0 0.0 0.0 0.000000
2014-12-14 0.0 0.0 0.0 0.000000
If the channel outlet has a rectangular shape of width 2m, how tall should our boots be? Assume a 2% slope, and find a reasonable assumption for the roughness of asphalt.
Recall the Manning equation:
Where:
n is the manning roughness
A is cross sectional area of the flow
R hydraulic radius (area / wetted perimeter)
S is the channel slope
w_channel = 1.5 # m
S = 0.005 # channel slope
n_factor = 0.017 # rough asphalt
def calc_Q(d, w, S, n):
"""
Calculate flow from the Manning equation.
"""
A = d * w
wp = w + 2 * d # wetted perimeter
R = A / wp
return (1/n) * A * R**(2/3) * S**(1/2)
def solve_depth(w, n_factor, Q, S):
"""
Given a flow, a roughness factor, a channel slope, and a channel width,
calculate flow depth.
"""
e = 1 / 100 # solve within 1%
d = 0
Q_est = 0
n = 0
while (abs(Q_est - Q) > e) & (n < 1000):
Q_est = calc_Q(d, w, S, n_factor)
# print(Q, Q_est, abs(Q_est - Q))
d += 0.001
n += 1
# print('solved in {} iterations'.format(n))
return d
# For each timestep, we want to solve for the depth of water at our outlet
sample_df['flow_depth_m'] = sample_df['runoff_cms'].apply(lambda x: solve_depth(w_channel, n_factor, x, S))
plt.plot(sample_df.index, sample_df['flow_depth_m'])
plt.ylabel('Flow depth [m]')
Text(0, 0.5, 'Flow depth [m]')
Not only are our feet wet, but if we happen to be there the peak it’s potentially dangerous. As little as 10-15cm of water moving fast enough can sweep you off your feet.

More Complex Implementation: Spatial Data¶
As discussed in class, precipitation takes time to travel from where it fell to the basin outlet. Next we will estimate the runoff response in a real catchment, just upstream from the parking lot example in the FitzSimmons Creek basin.
Step 1: Instantiate a grid from a DEM raster¶
Some sample data is already included, but for extra data, see the USGS hydrosheds project.
# grid = Grid.from_raster('data/n45w125_con_grid/n45w125_con/n45w125_con', data_name='dem')
grid = Grid.from_ascii(path='../../data/notebook_5_data/n49w1235_con_grid.asc',
data_name='dem')
# reset the nodata from -32768 so it doesn't throw off the
# DEM plot
grid.nodata = 0
# store the extents of the map
map_extents = grid.extent
min_x, max_x, min_y, max_y = map_extents
Plot the DEM¶
NOTE: The cell below may take up to 30 seconds to load. Please be patient, it is thinking really hard.
The code below will plot the Digital Elevation Model (DEM).
Do you recognize any features of the terrain? Can you locate where it is?
Hover over the map (or touch if using a touchscreen) to see the coordinates in decimal degree units.
What does the precision of the coordinates represent?
i.e. what does 5 decimal places in decimal degrees equate to in kilometers?
You can interact with the plot by using the tools on the left (in vertical order from top to bottom):
pan: move around the map
box zoom: draw a square to zoom in on
wheel zoom: use the mousewheel (or pinch gesture on a touchscreen) to zoom in
box zoom: draw a square to zoom in on
tap: not yet implemented (but you can see the coordinates)
refresh: reset the map
hover: see the coordinates when hovering over the map with a mouse or pointer
# set bokeh plot tools
tools = "pan,wheel_zoom,box_zoom,reset,tap"
# show the precision of the decimal coordinates
# in the plot to 5 decimal places
TOOLTIPS = [
("(x,y)", "($x{1.11111}, $y{1.11111})"),
]
# create a figure, setting the x and y ranges to the appropriate data bounds
p1 = figure(title="DEM of the Lower Mainland of BC. Hover to get coordintes.", plot_width=600, plot_height=int(400),
x_range = map_extents[:2], y_range = map_extents[2:],
tools=tools, tooltips=TOOLTIPS)
# map elevation to a colour spectrum
color_mapper = LinearColorMapper(palette="Magma256", low=-200, high=2400)
# np.flipud flips the image data on a vertical axis
adjusted_img = [np.flipud(grid.dem)]
p1.image(image=adjusted_img,
x=[min_x], # lower left x coord
y=[min_y], # lower left y coord
dw=[max_x-min_x], # *data space* width of image
dh=[max_y-min_y], # *data space* height of image
color_mapper=color_mapper
)
color_bar = ColorBar(color_mapper=color_mapper, #ticker=Ticker(),
label_standoff=12, border_line_color=None, location=(0,0))
p1.add_layout(color_bar, 'right')
show(p1)
# -123.15512, 49.41293
#-123.14657, 49.41080
-123.14350, 49.40251